: NMDAR10.mod
: Jose Guzman, sjm.guzman@gmail.com
: Last change: Sat Oct 26 22:36:00 CEST 2013

TITLE  ten state kinetic NMDAR model 

COMMENT
-------------------------------------------------------------------------

Based on the 10 state model of NMDA-gated channel reported in
Kampa BM, Clements J, Jonas P, Stuart GJ (2004)
Kinetics of Mg2+ unblock of NMDA receptors: 
implications for spike-timing dependent synaptic plasticity. 
J Physiol 556:337-45

The state transition diagram:

    R-------A2R-------O
     \       |\        \
      \      | \        \
       \     |  \        \
        RMg-----A2RMg-----OMg
             |   |
             Df  |
             |\  |
             | \ |
             Ds \|
              \  DfMg
               \ |
                \|
                 DsMg

Note that agonist and Mg concentration are folded into rate constants,
instead of being treated as reactants themselves.

The code draws directly from Bjoern Kampa's NMDA_Mg.mod but
1.  is driven by events
2.  explicitly represents transmitter concentration by a
rectangular pulse with user-specifiable duration and amplitude.

-------------------------------------------------------------------------
ENDCOMMENT

NEURON {
    POINT_PROCESS NMDAR10 
    : states
    RANGE R, A2R, O, RMg, A2RMg, OMg, Df, DfMg, Ds, DsMg

    : parameter rates
    GLOBAL rb, rmb, rmu, rbMg, rmc1b, rmc1u, rmc2b, rmc2u
    GLOBAL Rb, Ru, Rd1, Rr1, Rd2, Rr2, Ro, Rc, Rmb, Rmu
    GLOBAL RbMg, RuMg, Rd1Mg, Rr1Mg, Rd2Mg, Rr2Mg, RoMg, RcMg
    GLOBAL Rmd1b, Rmd1u, Rmd2b, Rmd2u, rmd1b, rmd1u, rmd2b, rmd2u
    GLOBAL Rmc1b, Rmc1u, Rmc2b, Rmc2u

    : parameters of voltage-dependent Mg block
    GLOBAL valence, memb_fraction

    THREADSAFE : assigned GLOBALS will be per thread

    RANGE g, gMax, Erev, mg
    RANGE dur, nt
    NONSPECIFIC_CURRENT i
}

UNITS {
    (nA)   = (nanoamp)
    (mV)   = (millivolt)
    (pS)   = (picosiemens)
    (umho) = (micromho)
    (mM)   = (milli/liter)
    (uM)   = (micro/liter)
}

PARAMETER {
    Erev = 5   (mV)                 : reversal potential
    gMax = 500 (pS)   <0.0,1e9>     : maximal conductance
    mg   = 1   (mM)   <0.0,1e9>     : extracellular Mg concentration
    dur  = 1   (ms)   <0.0,1e9>     : duration of Glu in the cleft

    : parameters of voltage-dependent Mg block 
    valence = -2
    memb_fraction = 0.8

    : [Glu] dependent rates 
    Rb = 10e-3 (/uM /ms)    : receptor binding rate, R-> A2R
    RbMg = 10e-3 (/uM /ms)  : receptor binding rate in Mg, RMg -> A2RMg 

    : Mg independent rates 
    Ru  = 5.6e-6  (/ms)     : receptor unbinding,  A2R-> R
    Ro  = 10e-3   (/ms)     : receptor opening, A2R -> O
    Rc  = 273e-3  (/ms)     : receptor closing, O -> A2R 
    Rd1 = 2.2e-3  (/ms)     : fast desensitization, A2R -> Df
    Rr1 = 1.6e-3  (/ms)     : fast re-sensitization, Df -> A2R
    Rd2 = 0.43e-3 (/ms)     : slow desensitization, Df -> Ds
    Rr2 = 0.5e-3  (/ms)     : slow re-sensitization, Ds-> Df

    : Mg and voltage dependent rates
    Rmb = 0.05e-3 (/uM /ms)     : Mg binding open,  O -> OMg
    Rmu = 12800e-3  (/ms)       : Mg unbinding open, OMg -> O
    Rmc1b = 0.00005e-3 (/uM /ms): Mg binding close 1, R -> RMg
    Rmc1u = 2.438312e-3 (/ms)   : Mg unbinding close 1, RMg -> R
    Rmc2b = 0.00005e-3 (/uM /ms): Mg binding close 2, A2R -> A2RMg
    Rmc2u = 5.041915e-3(/ms)    : Mg unbinding close 2, A2RMg -> A2R
    Rmd1b = 0.00005e-3  (/uM /ms) : Mg binding to desensitization fast, Df -> DfMg
    Rmd1u = 2.98874e-3  (/ms)  : Mg unbinding to desensitization fast, DfMg -> Df
    Rmd2b = 0.00005e-3  (/uM /ms)  : Mg binding to desensitization slow, Ds -> DsMg
    Rmd2u = 2.953408e-3  (/ms)  : Mg unbinding to desensitization slow, DsMg -> Ds

    : Mg depedent rates
    RuMg  = 17.1e-3 (/ms)  : receptor unbinding with Mg, A2RMg -> RMg
    RoMg  = 10e-3   (/ms)  : opening with Mg, A2RMg -> OMg
    RcMg  = 548e-3  (/ms)  : closing with Mg, OMg -> A2RMg
    Rd1Mg = 2.1e-3  (/ms)  : fast desensitisation with Mg, A2RMg -> DfMg
    Rr1Mg = 0.87e-3 (/ms)  : fast re-sensitisation with Mg, DfMg -> A2RMg
    Rd2Mg = 0.26e-3 (/ms)  : slow desensitisation with Mg, DfMg -> DsMg
    Rr2Mg = 0.42e-3 (/ms)  : slow re-sensitisation with Mg, DsMg -> DfMg

}

ASSIGNED {
    v  (mV)  : postsynaptic voltage
    i  (nA)  : current, where i = g*(v-Erev)
    g  (pS)  : synaptic conductance
    nt (mM)  : glutamate concentration in the cleft
    onset (1): 1 if transmitter is in the cleft

    : calculations
    rb    (/ms)  : binding, [glu] dependent
    rmb   (/ms)  : blocking V and [Mg] dependent
    rmu   (/ms)  : unblocking V and [Mg] dependent
    rbMg  (/ms)  : binding, [glu] dependent
    rmc1b (/ms)  : blocking V and [Mg] dependent
    rmc1u (/ms)  : unblocking V and [Mg] dependent
    rmc2b (/ms)  : blocking V and [Mg] dependent
    rmc2u (/ms)  : unblocking V and [Mg] dependent
    rmd1b (/ms)  : blocking V and [Mg] dependent
    rmd1u (/ms)  : unblocking V and [Mg] dependent
    rmd2b (/ms)  : blocking V and [Mg] dependent
    rmd2u (/ms)  : unblocking V and [Mg] dependent

}

STATE {
    : Channel states (all fractions)
    R       : unbound
    A2R     : closed
    O       : open
    RMg     : unbound, with Mg
    A2RMg   : closed,  with Mg
    OMg     : open, with Mg
    Df      : desensitisized 1
    Ds      : desensitisized 2
    DfMg    : desensitisized 1, with Mg
    DsMg    : desensitisized 2, with Mg
}

INITIAL {
    R     = 1       : all in unbounded state
    nt    = 0 (mM)  : no transmitter in the cleft
    onset = 0       : 1 when transmitter is in cleft
}

BREAKPOINT {
    SOLVE kstates METHOD sparse
    g = gMax * O
    i = g * (v-Erev) * (1e-6) : current in nA
}

KINETIC kstates {

    rb    = Rb * (1e3) * nt 
    rbMg  = RbMg * (1e3) * nt 
    rmb   = Rmb * mg * (1e3) * exp((v-40) * valence * memb_fraction /25 (mV))
    rmu   = Rmu * exp((-1)*(v-40) * valence * (1-memb_fraction) /25 (mV))
    rmc1b = Rmc1b * mg * (1e3) * exp((v-40) * valence * memb_fraction /25 (mV))
    rmc1u = Rmc1u * exp((-1)*(v-40) * valence * (1-memb_fraction) /25 (mV))
    rmc2b = Rmc2b * mg * (1e3) * exp((v-40) * valence * memb_fraction /25 (mV))
    rmc2u = Rmc2u * exp((-1)*(v-40) * valence * (1-memb_fraction) /25 (mV))
    rmd1b = Rmd1b * mg * (1e3) * exp((v-40) * valence * memb_fraction /25 (mV))
    rmd1u = Rmd1u * exp((-1)*(v-40) * valence * (1-memb_fraction) /25 (mV))
    rmd2b = Rmd2b * mg * (1e3) * exp((v-40) * valence * memb_fraction /25 (mV))
    rmd2u = Rmd2u * exp((-1)*(v-40) * valence * (1-memb_fraction) /25 (mV))

    : kinetic reactions
    ~ R     <-> A2R     (rb, Ru)
    ~ A2R   <-> O       (Ro, Rc)
    ~ A2R   <-> Df      (Rd1, Rr1)
    ~ Df    <-> Ds      (Rd2, Rr2)
    ~ O     <-> OMg     (rmb, rmu)
    ~ RMg   <-> A2RMg   (rbMg, RuMg)
    ~ A2RMg <-> OMg     (RoMg, RcMg)
    ~ A2RMg <-> DfMg    (Rd1Mg, Rr1Mg)
    ~ DfMg  <-> DsMg    (Rd2Mg, Rr2Mg)
    ~ R     <-> RMg     (rmc1b, rmc1u)
    ~ A2R   <-> A2RMg   (rmc2b, rmc2u)
    ~ Df    <-> DfMg    (rmd1b, rmd1u)
    ~ Ds    <-> DsMg    (rmd2b, rmd2u)

    CONSERVE R + A2R + O + RMg + A2RMg + OMg + Df + DfMg + Ds + DsMg = 1
}

NET_RECEIVE (weight (mM)) {
    if (flag == 0) {
        if (onset == 0) {: receptor not active
            onset = 1
            nt = weight 
            net_send(dur, 1) : send flag==1 at dur
        }
        else { 
            net_move(t+dur) :receptor is active, then extend its duration
        }
    }
    else if (flag == 1){ :internal event (1) finish transmitter
        onset = 0
        nt = 0
    }
}
